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Abstract 

We have calculated the neutron electric dipole moment (NEDM) in the presence of the CP 
violating 9 term in lattice QCD with 2-flavor dynamical clover quarks, using the external electric 
field method. Accumulating a large number of statistics by the averages over 16 different source 
points and over forward and backward nucleon propagators, we have obtained non-zero signals of 
neutron and proton EDM beyond one standard deviation at each quark mass in full QCD. We 
have investigated the quark mass dependence of nucleon EDM in full QCD, and have found that 
nucleon EDM in full QCD does not decrease toward the chiral limit, as opposed to the theoretical 
expectation. We briefly discuss possible reasons for this behavior. 

PACS numbers: 11.30.Er, 11.30.Rd, 12.39.Fe, 12.38. Gc 
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I. INTRODUCTION 



A requirement for renormalizability in QCD allows a CP violating term with a free 
parameter 9qqd, thus called the 9 term, 

,2 



£<5CP — #QCD^|— oGfiuG^, (1) 



where G^ u and represent a gluon field strength and its dual, and g is a coupling constant 
of QCD. If the CP invariance is preserved in the strong interaction in Nature, #qcd must be 
zero in QCD. Indeed a current experimental bounds of the neutron electric dipole moment 
(NEDM) [l], 

d N < 6.3 x l(T 13 e • fm, (2) 

and the crude theoretical estimates , cLn = 0(10 2 )#qcd, gives the constraint that 

6*q CD < 0(1O~ 10 ). #qcd must be very small or even zero. 

In the presence of the electroweak interaction, however, the above conclusion on #qcd 
is modified. In the Weinberg-Salam(WS) model of the electroweak interaction, the quark 
mass term generated through Yukawa couplings by the spontaneous electroweak symmetry 
breaking is given by 

£ m = q l L (M C KM)ijqR + h.c, (3) 

with a quark mass matrix Mckm and left- and right-handed quark fields Qr^l- In order to 
transform the quark mass matrix to a real and diagonal form, U(l) as well as SU(Nf) chiral 
rotations are necessary, since Mckm is non-Hermitian in general. Through Adler-Bell- Jackiw 
anomaly, this U(l) chiral transformation shifts #qcd in eq.(pQ) to 

= Oqcb + arg det M C km- (4) 

Therefore the experimental bound (T5]) on leads to 6 < O(10~ 10 ): A subtle cancella- 
tion between a parameter in QCD (#qcd) and a phase of mass matrix in the WS model 
(arg det Mckm) should be fulfilled to keep the CP invariance in the strong interaction. This 
cancellation seems unnatural and thus a new mechanism must exist to explain the smallness 
of 9 ("strong CP problem"). 

A simplest solution to the strong CP problem is that one of quarks is massless, so that 
we can set 9 = by the chiral rotation of this quark without introducing complex phases to 
other quark masses. Detailed analyses by CHPT 5f| or lattice QCD [J, however, strongly 



indicate that up quark has a finite mass (1.5 < rtiu < 3.0 MeV Qj), therefore this solution is 
excluded. In the Peccei-Quinn (PQ) mechanism jg], 9 is promoted to a scalar field associated 
with a new symmetry, called PQ symmetry, which is slightly broken by the anomaly. The 
spontaneous PQ symmetry breaking automatically pick up an unique vacuum, in which 9 
vanishes. As a consequence of the symmetry breaking, a very light new scalar particle, called 
axion must exist in this mechanism. So far the axion, which is also one of the candidates 
for the cold dark matter, has not been observed yet, and both cosmological observations 
and accelerator experiments set a very narrow allowed region of the axion mass such that 
10" 6 < m a < 3 x 10" 3 eV [9]. 

An aim of our investigation in this paper is NOT to propose a new solution to the strong 
CP problem, but is to establish a reliable way of calculating the NEDM in lattice QCD. For 
small 9, the NEDM is proportional to 9 as d^ = dffd + 0(9 3 ). Various model calculations 
lead to different estimates, ranging that \d ( i ] \ = C(10~ 2 ~ KT 3 ) e-fm, and even its 
sign has not been determined yet. The lattice QCD has a potential to calculate d$ non- 
perturbatively and in a model independent way. Once a method of calculating the NEDM 
is established in lattice QCD for the case of a particular CP violating term given in (CO), it 



may become possible to extend NEDM ca 
such as the chromoelectric dipole moment 



culations to the case of new CP violating terms 



4( generated in SUSY models at high energy 10 1. 

Ill- 



Lattice calculations of d$ have remained to be notoriously difficult for a long time 
In our previous papers [13I ] we have investigated two methods of calculating d$ in 



quenched lattice QCD. In the first method [121 ]. we have calculated the NEDM form factor 
at non-zero momentum transfer, which becomes the NEDM in the limit of zero momentum 
transfer. In the second method, the NEDM has been directly extracted from the energy shift 
in the presence of the external electric field. In both cases, we have successfully obtained 
non-zero signals of EDM for neutron and proton. In this paper we extend the calculation of 
the NEDM to the case of 2-flavor full QCD, using configurations generated by the CP-PACS 
collaboration [ijj]. We employ the second method, the direct calculation of the NEDM with 
the external electric field, since no extrapolation of momentum transfer is needed, as opposed 
to the form factor method. This method, however, requires a large number of configurations 
to reduce statistical errors. Performing calculations at four different quark masses, we discuss 
a quark mass dependence of the NEDM in full QCD. In particular it is interesting to see 
whether the NEDM vanishes or not at zero quark mass as theoretically predicted. There 
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exists a previous study of the NEDM in lattice QCD with 2-flavor dynamical domain- 
wall quarks using the form factor method Unfortunately the signal of the NEDM is 
consistent with zero within a large error, so that only an upper bound is obtained for a value 
of 

The organization of this paper is as follows. In SecJTTl we give a definition of NEDM in 
the presence of uniform and static electric field. We then consider how this external electric 
field is introduced on the lattice, and discuss a violation of uniformity of the electric field 
due to boundaries of the finite lattice. We also explain our method of extracting the NEDM 
from the nucleon propagator. In Sec JIII|. simulation parameters for dynamical configurations 
are briefly summarized. Our main results are given in Sec lIVl Summary and discussions of 
this paper are presented in SecJVj In this paper we set a = 1 unless necessary. 

II. METHOD OF THE NEDM CALCULATION 

A. Definition of the NEDM 

— * 

In the presence of the constant and uniform electric field E, a change of energy for the 
nucleon state due to the 9 term is denoted as 

A£ CP = d% ] 6S ■ E + O{E 3 0, EO 3 ) (5) 

with the nucleon spin vector S. This leads to the following extraction of djy for E = 
(0,0, E z ): 

£ e + (E z ) - E 6 _{E Z ) = d$0E z + <D{E 3 Z 6, E z 6 3 ) (6) 

where £±(E Z ) is an energy of the nucleon state with S z = ±1/2 in the presence of the electric 
field E = (0, 0, E z ). Hereafter we simply denote as d^. 

B. Introduction of the electric field on the lattice and a boundary effect 

On the lattice, an external electric field E k is introduced into link variables via a replace- 
ment in the Wilson-Dirac operator with 

U k (x) — e e ^U k {x) = U 9 k (E,x), Ut(x) — er^E^x) = U q k (E,x), (7) 
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where e q is an electric charge of a quark flavor q. Note here that a complex factor i does not 
appears in the exponent since E k is defined in the "Minkowski" space while t is the time 
coordinate of the "Euclidean" lattice, and therefore U%(E,x) is no more unitary 1 . In our 
calculation we have made this replacement of the link variables only for the valence quark 
while gauge configurations have been already generated without this replacement in 
the quark determinant. This "approximation" is equivalent to ignoring the "disconnected" 
contribution from the electric field. It is noted that, at the first order of the electric field, 
this contribution vanishes for the 3 flavor QCD with m u = = m s , since the disconnected 
contribution does note depend on the quark flavor in the flavor SU(3) limit, and therefore 
is proportional to e u + Cd + e s = after summing over 3 flavors. 



As discussed in Ref . (13| , the introduction of the " Minkowski" electric filed destroys the 
periodic boundary condition of link variable U%(E, x), so that translational invariance of the 
electric field is violated at the temporal boundary. Indeed an effective electric fields, defined 
by E k (t E ) = {\nU%(E,t E + l) -\nU%(E,t E )}/e q with U k (t E ) = 1, becomes 

E k at t E = 1,2, •••T- 1 
E k {t E ) = , (8) 

— (T — l)E k at t E = T 

where t E runs from 1 to T, so that t E = T + 1 is equal to t E — 1. A strong anti-electric field 
is generated between t E — 1 and t E = T in order to cancel the constant electric field E k , 
so that Y^t =i Ek{t E ) vanishes. We denote the place where the large gap exist as tca P , and 
^Gap = in the above case. Since the definition of the NEDM in ([61) requires the constant 
electric field, the NEDM should only be measured far away form the boundary to avoid an 
effect of the anti-electric field. In this work we keep distances between source/sink points 
and the boundary as large as possible. In particular we fix |tcap — ^src + 1| —T/2, where t src 
is the time coordinate of the source point. 



In the previous study 13] we found that the good sampling of the topological charge is 
important for obtaining the reliable signal of NEDM. At least an order of a few thousands 
configurations was needed to realize the symmetric and Gaussian distribution in quenched 
QCD. On the other hand, a number of full QCD configurations in this study is limited to 



1 With the electric field in the Minkovski space, however, the energy difference [5] becomes real in the 
lattice calculation, so that it can be extracted from the ratio nucleon propagators between spin-up and 



spin-down 11 1 
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700 ~ 750 at each quark mass. We therefore use several source points for each configuration 
to increase statistics. In accordance with the change of the source point, we shift the 
boundary point, keeping the distance between the source point and the boundary as large 
as possible, to avoid the boundary effect mentioned above. 



C. Spinor structure of the nucleon propagator and EDM 

In the presence of the electric field E and the CP violating 9 term, the explicit form of 
the nucleon propagator becomes 



(N a Np)g(E,t) = Z N (E 2 )[(l + A N (E 2 )9a ■ E 



x exp(-S N (E 2 )t-^-a-Et)} +••■, (9) 

2 Jo/3 

with an overall amplitude Z^(E 2 ), a coefficient A^(E 2 ) and an energy Sn(E 2 ) 2 . The ellipse 
represents O(E 3 0, 9 2 ) terms and contributions from excited states. Here (O)e represents a 
vacuum expectation value in reweighting method as 

Ze(O} = Z e=0 (Oe WQ ) e=0 (10) 

with topological charge Q, which is evaluated by the cooling method in our calculation 



13|- 
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In order to extract from the nucleon propagator in the above equation, we consider 
the following ratio between different spinor components; 

Rf nc {E, t; 9) B% aive (E = 0, t; 9 = 0) 



#naive(£ = Qj t] 9) R™ VC (E, t\ 8 = 0) 

1 + 9A N (E 2 )E 



1-9A N (E 2 )E 

where 



exp[-d N 6Et], (11) 



m t; 9) = V V A 'h: ,U, - AV ; = l + A ; 6 } E T e M -d N 9Et] + ■ ■ ■ , (12) 
3 K J (N 2 N 2 )e((0,0,E),t) 1-A N 9(E 2 )E PL N J ' K J 

Three additional i?3 aive 's in Rs{E, t; 9) are introduced to remove contamination coming from 
9 = and/or E = terms due to the insufficient statistics. In addition, to remove fictitious 



2 Note that because of the acceleration for a charged particle the energy of proton increases as a time 
increases [l6| . However this contribution is canceled out in the following ratio. 
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E 2n 9 contribution, we construct the ratio 

R^(EfO) R ^ E ^ 6 ^ Kr°{E,t;B) i^-fl,*; g = 0) 
3 1 ' ' ' R 3 (-E,t;9) R^ ve (-E,t;9) R™ ve (E, t; 9 = 0) 

with R 3 (E,t;9) in eg. (1111) . For later use we define the "effective djv" as 

'R c 3 OTI (E,t- 1-9 



2d N 9E = In 



(14) 



Rf TT (E,t; 9) 

which is validated in large t where the nucleon asymptotic state dominates, 

III. SIMULATION PARAMETERS 

In our calculation, we employ 2-flavor dynamical QCD configurations, generated by CP- 
PACS collaboration 14| with the RG-improved (Iwasaki) gauge action and the clover quark 
action on 24 3 x 48 lattice at /3 = 2.1. The corresponding lattice spacing, determined by the 
rho meson mass m p = 768.4 MeV, is a -1 ~ 1.8 GeV (a ~ 0.11 fm). 

The valence quark mass is chosen to be same with the sea quark mass. The pion mass 
m PS becomes 1.13, 0.93, 0.76 and 0.53 GeV at K sca = 0.1357, 0.1367, 0.1374 and 0.1382, 
respectively. Table fl] lists lattice parameters used in our simulation. Throughout this paper 
statistical errors are estimated by the jack-knife method, whose bin size is 5 configurations, 
equivalent to 25 HMC trajectories. We employ a local sink and a smeared source for all 
three quark propagators in the nucleon 2-pt function, with the exponential smearing, f(r) = 
Ae~ Br . Parameters (A, B) depend on the quark mass as shown in Table [B 

We take 16 different source points separated by 3 lattice units in temporal direction and 
maximally separated in spatial directions. Averaging these source points, a total number of 
statistics is more than 10,000. 

Throughout our study, the value of electric field and 9 are fixed to E = 0.004 and 
9 = 0.025, which are the most suitable choice to reduce statistical errors [l^ |. 
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IV. RESULTS 



A. Topological charge distribution and nucleon mass 

We measure the topological charge of each configuration after 50 cooling steps, using the 
O (a 2 ) improved definition of the topological charge density given by 



with a = 5/3 and a\ = -1/12 JlTj, where and L 1 ^ 2 are 1 x 1 and 1x2 Wilson 
loops, respectively. We show the time history of the topological charge in FigfT] and its 
histogram in Fig j2] at four quark masses. The distributions of the topological charge is more 
or less gaussian at all quark masses. As the sea quark mass decreases, the width of the 
gaussian distribution becomes narrower in accordance with the theoretical expectation. 

The effective mass of the nucleon at four quark masses is plotted in Figj3l The average 
over 16 sources gives enough statistics to produce a clear plateau at all cases. We see that 
the plateau for the nucleon state starts at t — t BIC + 1 = 6. The global fit of the propagator 
from t — t src + 1 = 9 to 15 gives a value of the nucleon mass with a very small error, as shown 
in Table [H 

B. Signal of EDM 

In the previous study we found that the influence due to the gap of the electric field 
on the EDM signal disappears at |tcap ~ t BIC + 1| > 5 at a = 0.1 fm in quenched QCD 
jl^ . Since we fix |tcap — ^src + 1| = T/2 = 24 in this study, the above condition is well 
satisfied. In this set up, the backward propagation of the nucleon becomes identical to 
the forward propagation in the infinite statistics, so that we can take an average over both 
to increase statistics. In FigJH we compare the time dependence of Rs(E,t;9) between 
forward and backward propagation for neutron and proton. We have found that two results 
marginally agree with each other at t — t src + 1 < 10. Here we consider that the difference at 
t — t BIC + 1 > 10 is caused by the statistical noises, and thus the signal of EDM is obtained 
only at t — t src + 1 < 10. Therefore a fitting range should be chosen as 6 < t — t src + 1 < 10. 
Since T — (t — t Gap ) = T/2 — 9 = 15, the distance between the gap and the end-point of the 
fitting ranged at t — t STC + 1 = 10, is much larger than 5, the gap of the electric field does 



improved 



1 

32vr 2 




(15) 



S 



not affect the EDM signal. 

In Figj5JlHl the time dependence of Rs(E, t; 9) with E = ±0.004 is plotted for the neutron 
and proton at each quark mass. At all quark masses, the signal for non-zero EDM can be 
seen at 6 < t — t STC + 1 < 10, and the signal changes its sign as E does. To estimate the size 
of nucleon EDM, we consider an effective mass of Rl°"(E,t; 9) defined in eq. fll4p . and plot 
it in Figj9lfT2l Although errors and fluctuations are large, non-zero signals for both proton 
and neutron EDM have been observed in full QCD simulations for the first time. By fitting 
R^°"(E, t; 9) with an exponential function in eq. (fT3|) at 6 < t — t STC + 1 < 10, we obtain the 
value of d$ , which is given in Table fll] Signs of EDM are opposite between proton and 



neutron. This agrees with the quenched result [13J and with the CHPT prediction 



a. 



C. Mass dependence of the nucleon EDM 

In FigJTBl for neutron and proton are plotted as a function of pion mass squared, mp S , 
together with the quenched results. Unfortunately, because of large statistical errors in full 
QCD results, it is hard to observe a difference from the quenched results. Compared with 
the model calculation J2J, the central value is 10 times larger although its errors are large. 

Mass dependence of d^ in full QCD does not show the expected decrease toward the 
chiral limit. There are several possible explanations. Firstly, large statistical errors might 
hide the actual decrease of d^ toward the chiral limit. Secondly, the quark mass in this full 
QCD simulation is still too heavy to see the decrease. Thirdly, d^ does not vanish in the 
chiral limit due to the lattice artifact that the chiral symmetry is explicitly broken in the 
Wilson-type quark action. Indeed the topological susceptibility, which is also theoretically 
expected to vanish in the chiral limit, does not show the decrease in full QCD configurations 
of this paper [3]. In FigJHJ d^ is plotted as a function of the topological susceptibility 
instead of the pion mass squared. Contrary to the case of the pion mass squared, the change 
of the topological susceptibility is too little to observe a possible decrease of djy with large 
statistical errors. This suggests that the combination of the 2nd and 3rd possibilities is a 
main reason for the mass dependence of d^ in Fig. [131 

We also present the results of the mass dependence for the nucleon CP-odd phase factor 

nn 

in full QCD calculation. As discussed in |12l . I X 31 ] . the next-leading term in the nucleon 
spinor structure contains additional phase factor f^, which arises from m^e 4 ^ 975 = m^il + 
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if N e l5 ) + o(9 2 ), as 

(N(p,t)N(p,0)Q) = l^| 2 e-^4^T5- (16) 

where p, Zn, En denotes the nucleon momentum, amplitude and energy, respectively. This 
factor should goes to zero toward the massless limit because of the same reason as NEDM 
does. In Fig. [151 fk m QCD is plotted as a function of rrip S , together with the quenched 
results. Although the statistical errors are not so large compared to the EDM results, it 
does not show the expected decrease toward the chiral limit as well as the EDM case. This 
observation suggests that it is unlikely that the correct chiral behavior of the EDM is hidden 
in its large statistical errors. 



V. SUMMARY AND DISCUSSION 



In this paper we present the evaluation of the NEDM in full QCD simulation using the 
external electric field method. After accumulating a huge number of statistics with multiple 
sources and an average over forward and backward propagation, we have obtained non-zero 
value of the nucleon EDM for the first time in full QCD. Statistical errors of the EDM in 
full QCD are still larger than in the previous quenched case. The mass dependence of the 
EDM, which is similar to the quenched one, does not show the expected behavior in full 
QCD that it vanishes towards the chiral limit. Besides the large statistical errors, there may 
be two main reasons. One is that sea quark masses used in this paper are still too heavy 
to see the expected decrease, the other is that the explicit chiral symmetry breaking of the 
Wilson type quark action spoils the expected chiral behavior. Indeed this chiral behavior of 
the EDM is very similar to that of the topological susceptibility. 

Since it now becomes possible to calculate the EDM in full QCD, we should proceed 
toward the precise evaluation of the NEDM. First of all, we should further decrease the 
sea quark mass to clearly observe the chiral behavior of EDM. As reported in the recent 
advanced work in PACS-CS collaboration 19|, 2 + 1 flavors dynamical configurations are 



being generated at the quark mass close to the physical point, and these configurations 
will be available soon. Secondly, we should decrease statistical errors of the EDM. For 
this purpose, it may be better to switch to the form factor method 12J, though the zero 
momentum transfer limit has to be taken in this case. From the previous work 13J and this 
work, the good chiral behavior of the quark action seems not so relevant to obtain the signal 
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of the EDM, and we are currently calculating the EDM form factor using the clover quark 



action 



201 ] . Finally the effect of disconnected loop diagram, ignored in our studies, should 



also be included in the calculation. A preliminary study shows that it is possible to include 



this effect with the current available resources 
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#configs. 

FIG. 1: Time histories of the topological charge at each sea quark mass. 
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FIG. 3: The nucleon effective mass as a function of time in lattice unit without electric field. The 
average over 16 source sets on each configuration is taken. 
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FIG. 4: Comparison between forward and backward propagation of R3 for neutron as a function of 
time at E = ±0.004 and = 0.025. The average over 16 source sets (t src = 3, 6, • ■ ■ ,48) is taken. 
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FIG. 5: -R3 as a function of time at K = 0.1357, after averaging 16 source sets (t STC = 3, 6, • • • , 48), 
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FIG. 11: Same as Fig. Mat K = 0.1374. 
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FIG. 12: Same as Fig. Mat K = 0.1382. 
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Tables 



TABLE I: Table for the simulation parameters in this work. 0{a) improved coefficient of the clover 
term csw is chosen as csw = 1-47. The column of (A, B) denotes the smeared source parameters, 
and K denotes the hopping parameter for the degenerate up and down quarks. 

Lattice size Physical volume (fm 3 ) cutoff a" 1 (GeV) K (A, B) mps/my mxa 
24 3 x 48 2.6 3 1.83 0.1357 (1.5,0.45) 0.81 1.1851(10) 

0.1367 (1.5,0.43) 0.76 1.0224(11) 
0.1374 (1.5,0.35) 0.69 0.8944(12) 
0.1382 (1.5,0.25) 0.58 0.7167(14) 



TABLE II: Results for neutron and proton EDM at each quark mass 



K mps (GeV) neutron EDM (e-fm) proton EDM (e-fm) 
0.1357 1.13 -0.014(11) 0.049(21) 

0.1367 0.93 -0.046(16) 0.019(29) 

0.1374 0.76 -0.031(16) 0.060(28) 

0.1382 0.53 -0.040(28) 0.072(49) 
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